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Abstract 

We consider several patchy particle models that have been proposed in literature and we inves- 
tigate their candidate crystal structures in a systematic way. We compare two different algorithms 
for predicting crystal structures: (i) an approach based on Monte Carlo simulations in the isobaric- 
isothermal ensemble and (ii) an optimization technique based on ideas of evolutionary algorithms. 
We show that the two methods are equally successful and provide consistent results on crystalline 
phases of patchy particle systems. 
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I. INTRODUCTION 



Over the last number of years considerable effort has been dedicated to predict the crys- 
talline phases for a wide variety of model systems. In the case of strongly interacting 
systems, such as atomic and molecular ones, much of the phase behavior is governed by the 
zero temperature case. In such situations, techniques which minimize the thermodynamic 
potential (e.g., genetic algorithms or Monte Carlo basin hopping simulations) have proven 
to be very useful in predicting the ground state structures^-. Motivated by these successful 
approaches, the aforementioned optimization procedures have been extended to soft mat- 
ter systems and turned out to be efficient and robust techniques, suitable for a reliable 
prediction of crystalline structures at zero temperature for a broad variety of systems^ - —. 
Optimization strategies search for particle arrangements that minimize the thermodynamic 
potential of a system at zero temperature and identify them as candidate equilibrium struc- 
tures. For finite temperatures this minimization criterion cannot be completely trusted as 
one cannot safetly neglect the entropic contributions to the thermodynamic potential. In- 
deed, it has been shown that, even though the ground state structures provide a good guess 
for the finite temperature candidate structures, crystal phases that are local minima of the 
thermodynamic potential can be thermodynamically stable at finite temperature*^. Opti- 
mization techniques which combine quasi-Newton local and global optimization steps can 
be successfully applied to soft matter systems thanks to the smoothness of the inter-particle 
interaction, which is needed to guarantee that the derivatives involved in the minimization 
procedures of the potential energy are continuous. 

In the case of hard (colloidal) particles, the situation is considerably different and much 
more difficult. The phase behavior in purely hard systems is completely governed by the 
entropic contribution to the free energy. Hence, when applying minimization techniques to 
hard systems in order to identify the stable ordered equilibrium structures, the question 
regarding what to minimize arises. Frequently in the past, the maximum packing fraction 
criterion has been used as it minimizes the Gibbs free energy at infinite pressure^ - — . For 
finite pressures, the entropic term cannot be neglected and hence this criterion is not fully 
reliable. For instance, for binary hard-sphere mixtures, crystalline structures which are not 
the best packed ones exist as stable phases in the phase diagram^. 

Recently, a method has been proposed to cope with hard-particle systems^. This ap- 
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proach is a statistical sampling method used as a search strategy for candidate crystals 
structures of given systems. It is based on simple Monte Carlo simulations of only a small 
number of particles: simulations are carried out at constant temperature and pressure in a 
simulation box, whose shape is free to fluctuate^ 1 ^. In Ref.— , Filion et al. demonstrate 
that the method successfully predicts both the infinite pressure as well as the finite pressure 
crystal phases for a variety of hard-core systems including binary hard-sphere mixtures, 
oblate hard sphero cylinders, hard asymmetric dumbbells and hard bowl-shape particles. 
Additionally, the method turned out to be successful for systems where particles interact 
via long range Coulomb interactions and Lennard- Jones interactions. Hence, the variable 
box shape approach also offers the possibility of predicting crystal phases at finite tempera- 
tures, as opposed to the zero temperature minimization techniques. However, we emphasize 
that these techniques determine only candidate crystal phases: subsequent full free energy 
calculations are still required to identify the stable phases and to draw the complete phase 
diagram^iI&I&IZ. 

The present paper is dedicated to the prediction of crystal structures of patchy particle 
systems. To be more specific, we study a variety of patchy models that have been pro- 
posed in literature^ and search for their crystal structures, using both the variable box 
shape simulation method as well as an optimization approach based on ideas of evolution- 
ary algorithms. Our motivation resides in the fact that patchy particles have become a 
class of promising colloidal particles that are able to self-assemble as building units of fu- 
ture materials*^— , with a host of wide-spread applications, ranging from photonic crystals 
to biomaterials. Thus, controlled synthesis and abundant production of colloidal particles 
carrying a specific (chemical or physical) pattern on their surfaces has become a hot topic 
in soft matter physics^ - — . A reliable prediction of the ordered equilibrium phases of such 
systems represents therefore an important element in designing larger, functional units with 
desired properties. 

While, for patchy particles, an evolutionary algorithm approach has recently demon- 
strated its power in successfully predicting candidate crystal phases that have shown to be 
stable at finite temperatures^^, in this contribution we apply for the first time the variable 
box shape technique to model systems with strong directional interactions. Hence, in the 
body of the paper we discuss candidate crystals obtained via variable box shape simulations 
for different models of patchy systems. Our results are compared -whenever possible - with 
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the available literature, in an effort to validate this method. For some of the investigated 
models, we directly compare results from the variable box shape simulation technique to 
those from the evolutionary algorithm method. We show that the two methods provide 
consistent results for the investigated patchy models. 

The paper is organized as follows. In section [III we introduce the selected models for 
patchy particles. In section IIIH we briefly describe the numerical methods, referring for 
more details to previous publications. Section [IV] is dedicated to the discussion of the 
ordered crystal structures for a selection of different patchy particle models. Finally, in 
section [V] we draw our conclusions and provide a comparison between the two algorithms 
when applied to predicting ordered phases of patchy systems. 

II. MODELS 

We consider several patchy particle models that have been proposed in literature during 
the past years. We distinguish two main categories: discontinuous, square-well type poten- 
tials and continuous, Lennard- Jones type interactions. For all the investigated models we 
focus on the single bond per patch regime. 

A. Orientational Square- Well Models 

Most of the patchy particle models proposed and used in literature are based on hard-core 
particles whose surfaces are decorated by a fixed number of bonding patches interacting via 
a square-well potential. Here, we consider two orientational square-well models: the "sticky 
spots" model^ and the Kern-Frenkel model^i. 

The sticky spots model consists of hard-spheres carrying a small number of attractive 
points arranged in a regular geometry on the particle surface. The pair potential between 
two particles is given by the sum of an isotropic hard-core repulsion of diameter a and a 
site-site attraction. Sites on different particles interact via a square-well potential of depth 
e and attraction range 5 = 0.119cr; this choice of S guarantees that each site is involved 
at most in one bond^. In Section IIVI we show results for the case of six sticky sites per 
particle, thus each particle can form up to six bonds. Consequently, the average energy per 
particle, e, can vary from (system of monomers) down to — 3e (fully bonded system). 
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In the Kern-Frenkel model, pairs of particles interact via a square- well potential of depth e 
and attraction range 5, modulated by a function that depends only on the relative orientation 
of the attractive patches located on the two interacting particles. This function is zero when 
the patches do not feel each other; in this case particles only experience the repulsion 
due to the hard-core of diameter a. Otherwise, when the patches can feel each other, the 
modulating function is equal to unity Two patches on different particles feel each other when 
(i) particles are separated by a distance smaller than a + 5, and (ii) the vectors connecting 
the center of the particle with the center of the patch form an angle less than a maximum 
angle, fl max , with the vector connecting the centers of the two particles. By appropriately 
choosing 9 max and 5, multiple bonds between patches can be avoided and the single bonding 
regime is guaranteed^. In Section IIV[ we show results for the case of particles with four 
patches arranged in a tetrahedral geometry. The patch-patch attraction range and the patch 
angular extension are chosen to be S = 0.24cr and # max = acos(0.92), respectively. The phase 
behavior of this particular system has been investigated in Refs.— 

B. Orientational Lennard-Jones Models 

More realistic models for patchy particles describe the directional pair interaction via 
continuous and smooth pair potentials, which are in general longer ranged than their orien- 
tational, square-well counterparts. We consider a patchy particle model first introduced in 
Ref.— . In the model, the repulsion between two particles is given by an isotropic Lennard- 
Jones repulsive core, while the directional patch-patch attraction is specified by a Lennard- 
Jones attraction of depth e, modulated via a Gaussian-shaped angular decay. Provided that 
the patches are sufficiently narrow, the single bonding condition is guaranteed. We choose 
the following interaction parameters: the cut off of the attractive tail is r c = 2.5<r and the 
width of the Gaussian modulation is u = 0.3 rad^ 1 ^. In Section IIVI we show results for 
the cases of four and six patches per particle arranged in a tetrahedral and an octahedral 
geometry, respectively. The phase diagrams for such systems have been investigated in 
RefsM 

Other geometrical arrangements of the patches within the orientational Lennard-Jones 
pair potential can be introduced via an additional geometrical parameter g, which controls 
the patch positions on the particle surface^. Here we focus on the case of four patches per 
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particle, with r c = 1.9a and u = 0.22 rad. The parameter g is the central angle between a 
patch chosen as the pole of the particle and any of the other three patches. Hence, g ~ 109° 
specifies a regular tetrahedral arrangement. By varying g, the arrangement of the four 
attractive patches ranges from a rather compressed to a rather elongated tetrahedron, so 
that it is possible to study the effect of the patch geometry on the crystal structures. 

III. METHODS 

To predict crystalline structures, we use two different methods: the variable box shape 
simulation method of Ref.— and an optimization approach based on ideas of evolutionary 
algorithms'^. 

The identification of the various ordered structures and the comparison of the results 
obtained from the two different algorithms are based both on visual inspection and on a 
bond-order parameter analysis^ 1 ^. 

A. Variable Box Shape Simulations 

The algorithm is based on Monte Carlo (MC) simulations in the isothermal-isobaric 
ensemble: the initial state point is chosen to be in the fluid phase and then, at constant 
temperature T, the pressure P is increased step by step until the system converges to a final 
crystalline structure. We treat the simulation box as a unit cell and we allow its shape to 
fluctuate in order to avoid any bias of the crystal structure. As a consequence of the chosen 
box sizes, we work with extremely small numbers of particles N. In this paper, the particle 
numbers in the simulation box range typically from 4 to 8, but we also run simulations with 
up to 12 particles in the unit cell. Each MC step consists on average of N trial particle 
moves and one trial volume change, where the corresponding acceptance rules are given by 
the Metropolis algorithm. A particle move is defined as both a displacement in each direction 
of a random quantity distributed uniformly between ±5r and a rotation around a random 
axis of a random angle distributed uniformly between ±56. A volume move is given by a trial 
change of a randomly chosen component of a randomly chosen vector of the simulation box 
by a random quantity uniformly distributed between ±5v. The chosen values for the trial 
changes are 5r = 0.05<r, 56 = 25r/a rad, and 5v = 0.02cr, but they are allowed to change 
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during the simulation runs in order to keep the acceptance rates of both types of trial moves 
between 30% and 40%. The size of the trial moves is changed according to the following 
rule: if the acceptance probability -calculated every 10 5 Monte Carlo steps- of a particle (or 
volume) move is smaller than 30% then Sr = 0.955r (or Sv = 0.955?;); if the acceptance rate 
is bigger than 40% then Sr = 1.055r (or 8v = 1.055?;). Upper and lower limits for the step 
size are fixed in order to prevent extreme fluctuations in the fluid regime, namely 5r m i n (or 
5?; m in) = 0.001 and 5r max (or Sv max ) =0.5. Typically, both Sr and Sv increase abruptly in 
the fluid phase, because the initial configuration is in a very dilute regime; on progressively 
increasing pressure both Sr and Sv equilibrate fast to a value around 0.05 ± 0.05cr (bigger 
step sizes for higher temperatures). 

During the runs, we prevent the box from undergoing extreme distortions by using a 
lattice reduction technique^ in order to avoid extremely time consuming energy calculations. 
Moreover, we impose a lower bound on the length of all the lattice vectors: each of the lattice 
vectors must be longer than 1.5a. In this way, we avoid extreme elongations of the simulation 
box, in which the particles tend to form columns so that they only interact with the periodic 
images in one of the lattice directions. 

Once the number of particles, the temperature and the initial pressure are chosen, we run 
several MC simulations in parallel, starting from different initial conditions. The pressure 
is increased step by step, using on average 100 pressure steps from the initial to the final 
pressure; for each pressure value we perform 10 6 -10 7 MC steps. We distinguish two ranges 
for the final pressure: (i) low pressure values, ranging from 0.01 to 10 (in units of e/cr 3 ) and 
(ii) high/intermediate pressure values, ranging from 10 to 200. Different temperatures in 
the range from 0.10 to 0.20 (in units of e/ ' Kb) are considered. For each state point, we check 
if convergence to a certain final structure occurs over the last part (about one third) of each 
MC run. 

We note here that this method assumes that states which are stable for large systems are 
at least metastable for small systems, a point we feel to be largely validated by the fact that 
the method works well for the large variety of systems it has been tested for—. We also note 
that the small system sizes aid us in exploring phase space in two ways: (i) the fluctuations 
in density at a fixed pressure are larger than for larger systems; as such, near coexistence, 
the system frequently crosses the fluid-solid phase boundary and has a high probability of 
finding the stable phase, (ii) the small systems allow for large rearrangements of particles, 
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and hence significant changes in the crystal structure which would not be possible for large 
systems. 

B. Evolutionary Algorithm 

In an effort to identify the ordered equilibrium structures formed by the patchy particles at 
vanishing temperature we use optimization techniques that are based on ideas of evolutionary 
algorithms^ (EA). Working at fixed particle number (per unit cell) and fixed pressure P, 
the Gibbs free energy G is optimized with respect to (i) the lattice vectors specifying the 
unit cell, (ii) the positions of the particles within the unit cell, and (iii) the orientations of 
these particles. At zero temperature, G is reduced to the enthalpy: G — > H = U + PV, 
U being the lattice sum and V the volume of the system. In this contribution we use 
a phenotype implementation^ of such an algorithm, combining global optimization steps 
with local ones, as specified in Ref.— for the two-dimensional case. For technical details 
about the generalization to the three-dimensional case, we refer the reader to Ref.—. In the 
optimization runs, up to 8 particles per unit cell are considered. A population of usually ten 
individuals (each of them corresponding to an ordered candidate structure) is iterated along 
an evolutionary path via the usual mating, mutation and local minimization operations 
performed on the individuals (for details cf Ref.—). Throughout the optimization runs, the 
parameters of all these individuals are recorded. Among those, the one with the lowest 
value for the Gibbs free energy is considered as the final solution (global minimum) for this 
particular run; in addition, further structurally different^ local minima on the G-landscape 
are recorded. At least three and up to ten independent optimization runs are carried out in 
parallel for a given state point in order to ensure consistency. 

IV. RESULTS 

In the following, we discuss candidate crystals obtained via the variable box shape simula- 
tion approach for all the patchy models described in Sec. [TT1 Whenever possible, we compare 
our results with the available literature, in an effort to validate the method for systems char- 
acterized by strong directional interactions. Moreover, for the orientational Lennard- Jones 
models, we directly compare results from the variable box shape simulation technique to 
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those from the evolutionary algorithm method. All the lattice structures shown in the fig- 
ures are MC output data, while for the visual representation of the EA output data we refer 
to Refs.->^. The comparison between results from the two methods is reported in Table IIHI 
and FigJl 

A. Orientational Square-Well Models 

We first consider a particular realization of the Kern-Frenkel model with four patches, 
whose phase diagram has been extensively studied in Refs.— In these papers, the fol- 
lowing stable phases have been studied: the Face-Centered-Cubic (FCC) structure at high 
densities, the Body-Center-Cubic (BCC) crystal at intermediate densities, and two, tetrahe- 
drally arranged, open structures, i.e. the Diamond-Cubic (DC) and the Diamond-Hexagonal 
(DH) crystals^ 1 ^. In our MC simulations, we observe almost all the previously predicted 
crystal phases, only instead of the BCC lattice we identify a Body-Centered- Tetragonal 
(BCT) phase. In addition, we observe two Hexagonal-Close-Packed (HCP) crystals with 
different bonding patterns. 

Representative parts of all the above mentioned structures are shown in Fig. [TJ Since in 
the model the number of bonds is well defined, the bond saturation is indicated in the figure 
via a color code. For each structure, the corresponding values for the average energy per 
particle, e, and the average number density, p, are listed in Tab. HI The table also reports 
the frequency of appearance, / (expressed in percentage), of the structures encountered in 
the MC simulations: out of a total of 160 parallel simulations at high/intermediate pressure 
values, each of them starting with different initial conditions, 90% converged to one of the 
listed close-packed lattices. The corresponding values of / for the open structures turns out 
to be significantly smaller: out of a total of 70 simulations at low pressure values, only 53% 
converged to one of the two open configurations. This difference is due to the competition 
of the latter structures either with gel-like states or with hybrids between the DC and the 
DH lattices. It has been shown that, in large systems, hybrids of DC and DH structures are 
predominant^. 

As shown in Tab. (TJ the DC and the DH crystals have the same e- and p-values; however, 
in our simulation runs, we observe the DC structure with a slightly higher frequency than 
the DH lattice. Both diamond crystals are fully bonded structures built up of six-fold non- 
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planar rings. The difference between the two four-coordinated particle arrangements can 
be clarified by inspecting the bonds between adjacent layers: as highlighted in Fig. [1] by 
the yellow circles, pairs of particles forming intra-layer bonds occur for the DC crystal in 
the staggered conformation and for the DH case in the eclipsed conformation^. At higher 
densities, the close-packed structure with the highest /-value is a fully bonded FCC crystal. 
This lattice can transform into another fully bonded, but more compact lattice, the BCT 
crystal, which can also be viewed as a face-centered crystal with a non-cubic unit cell. 
Finally, the best packed crystals found for this model are two structures of HCP type with 
the same e- and p-values, but different bonding patterns; also the /-values of the two HCP 
lattices are considerably different. 

Another patchy model of the orientational square-well type is the sticky spots model 
introduced in Ref.— . Here, we consider particles decorated with six patches. To the best of 
our knowledge, the crystal phases of this model have not been investigated yet. In Fig. [2j 
we show the unit cells of the candidate ordered structures. Since in the model the number of 
bonds is well defined, we make again use of a color code to indicate the bond saturation of 
each particle. In Tab. UU we report the corresponding e-, p-, and /-values for each structure. 

As an open structure we consistently find the obvious Simple-Cubic (SC) crystal where 
all bonds are saturated. Additionally, we find (with a considerably lower /-value) another 
fully bonded structure, whose density is still smaller than that of the SC lattice. Such an 
ordered structure is built up of parallel, connected planes, in which particles are arranged 
in a honeycomb (Hcl) geometry, i.e. as six-fold planar rings. As candidate high pressure 
structures, we find HCP, FCC and BCT crystals. The structures with the highest densities 
are a partially bonded HCP structure (see panel (a) of Fig. [2]) and two FCC lattices (see 
panels (b) and (c) of Fig. [2]), one of which is fully bonded, while the other one is only 
partially bonded. The most frequently occurring high pressure lattice is a partially bonded 
FCC crystal (see panel (d) of Fig. [2]), whose energy is higher than that of the other partially 
bonded FCC crystal and whose density is significantly smaller than that of the two FCC 
lattices mentioned above. Finally, we identify the BCT crystal (see panel (e) of Fig. EJ) as a 
fully bonded structure with a relatively high density. 

We note that the percentage of simulations that converged to one of the close-packed 
structures listed in Tab. [Ill adds up to 57%, while the corresponding total /-value for the 
open lattices is 50.6%. The /-value for the low pressure simulations is comparable to the 
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Kern-Frenkel case discussed above, indicating once more the competition between the open 
lattices and gel-like states. For the close-packed lattices, instead, the considerably lower 
success rate of the present model as compared to the Kern-Frenkel model is related to the 
abundance of FCC structures with e varying from the fully bonded case, i.e. e = — 3e (see 
panel (b) of Fig. [2]) to e = — 2e (see panel (d) of Fig. [2]). An example of a FCC structure 
with an intermediate e-value is shown in panel (c) of Fig. |2j 

B. Orientational Lennard- Jones Models 

Next, we consider the orientational Lennard- Jones model with four and six patches per 
particle. In our MC simulation, we succeed to identify all the structures reported in liter- 
ature^ 1 ^. For patchy particles carrying six patches on their surface, we find all and only 
the predicted phases, i.e. FCC, BCC and SC. The corresponding unit cells are depicted in 
Ref.—. In the four-patch case we identify the FCC and BCC phases reported in Ref.—. In 
addition, we observe the DC and the DH lattices. The corresponding particle arrangements 
of the observed structures are similar to the ones displayed in panels (d)-(f ) of Fig. [TJ It has 
been showr*2& that the stability of the diamond structure sensitively depends on the posi- 
tion of the potential minimum (i.e., on the optimal distance between two bonded particles) 
and consequently on the attractive interaction range. For the chosen potential model, both 
diamond crystals are not thermodynamically stable^. This is an evidence of the need for 
full free-energy calculations to investigate the stability of the candidate crystal structures 
found. 

We also consider a related orientational Lennard- Jones model with four patches arranged 
in different tetrahedral geometries on the particle surface^, specified by the geometrical 
parameter g. Here, we show results for two extreme patch arrangements, i.e. g = 90° 
(compressed tetrahedron) and g = 150° (elongated tetrahedron) (see left column of Fig. [3] 
for a schematic representation). We compare candidate structures proposed by the variable 
box shape MC simulation technique and lattices suggested by the evolutionary algorithm 
approach. The candidate structures found via MC simulations are shown in Fig. [3J the 
corresponding e-, p-, and /-values of the structures are listed in Tab. IIHI 

In the case of compressed tetrahedrons (g = 90°), we observe, at low pressure values, 
the formation of either honeycomb double layers (HcDl) or hexagonal double layers (HxDl). 
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The HcDl structure is characterized by a fully bonded, planar honeycomb lattice, where the 
bonds within the six-fold rings are formed by the patches located on the equatorial plane 
of the particles. The remaining patches of each of the six particles forming a ring are all 
oriented in the same direction, providing the intra-layer link: two oppositely oriented layers 
are connected via these intra-layer linkers, forming hereby double layers. Between the double 
layer, there is no attractive interaction, thus they can either be far from each other or they 
can be almost in contact (the density value reported in Tab. Ill II refers to the latter case). 
The HxDl structure can be viewed as a HcDl structure with an additional particle located in 
the center of the six-fold honeycomb ring. In this case, double layers do interact with each 
other since the central particles provide links between the double layers, leading thereby 
to a higher p-value as compared to the preceding case. The e-value for the HxDl lattice 
is considerably higher than that of the HcDl, since, in this particle arrangement, the six- 
fold rings are slightly distorted in order to appropriately accommodate the seventh particle 
in their center. At high/intermediate pressure values, the most frequently encountered 
structure is the FCC crystal. A slightly more compact structure that is observed is a HCP 
crystal with higher e- and p- values. We note that, for this particular model, the /-values 
are throughout higher than those reported for the square-well patchy models. To be more 
specific, the percentage of simulations that converged to any of the packed structures at 
high/intermediate pressure values amounts to 97%, while the percentage of low pressure 
runs that converged to any of the open structures at low pressure values is 76%. 

Among the structures obtained from the MC simulations, only two are identified as 
global minima of the Gibbs free energy at zero temperature by the evolutionary algorithm 
method: the HcDl, at low pressure values, and the FCC crystal, at high pressure values. 
The additional crystals obtained via the MC simulations, i.e. the HxDl and the HCP 
structures, are identified as low-lying local minima on the Gibbs free energy landscape by 
the evolutionary algorithm. In Tab. IHIl we report the corresponding e- and p-values of 
the structures identified by the evolutionary algorithm together with the information about 
the location of the corresponding minimum on the Gibbs free energy landscape at zero 
temperature. In the intermediate range of pressure, the evolutionary algorithm identifies 
two additional global minima corresponding to two different kinds of hexagonal double 
layers; they are not found via the MC simulation technique. In Tab. IHIl we refer to them 
as HxDl-I and HxDl-II. For the visualization of such hexagonal double layers, see Ref.— . 
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In the case of elongated tetrahedrons (g = 150°), we observe, at low pressure values, the 
formation of a tetrahedrally arranged structure (Ts) (see panel (h) of Fig. [3]) as well as the 
formation of a slightly more compact lattice built up of staggered and connected double layers 
with a hexagonal structure (HxSDl) (see panel (g) of Fig. [3]), which is characterized by a 
slightly lower energy than the Ts configuration. The bonding patterns of the two lattices are 
distinctively different. In the Ts lattice each patch on a particle is strongly interacting with 
a patch on a neighboring particle: in the part of the Ts structure shown in Fig. [31 the upper 
particle and one of the lower particles are oriented with the polar patch pointing upwards, 
while the other two particles at the bottom are oriented downwards. In the HxSDl lattice, 
instead, pairs of oppositely oriented layers are bound to each other via the polar patches 
of each of the seven particles forming the hexagonal tiling; however, among the remaining 
patches, only two of them (per particle) are strongly interacting with the corresponding 
patch on the other, oppositely oriented layer. At high and intermediate pressure values, we 
identify either FCC or HCP crystals: the FCC structure is the most probable one; on the 
other hand, the HCP crystal maximizes the density but has a significantly lower e-value in 
comparison to the FCC lattice. Once more, we note that, for this model, the frequencies 
of occurrence sum up to values significantly higher than those reported for the square-well 
type patchy models: 99% for the packed structures and 75% for the open lattices. 

As we compare the results from the two different methods, we note that all the global 
minima identified by the evolutionary algorithm at zero temperature, i.e. the HCP and the 
HxSDl structures, are also obtained from the MC simulations at low (but finite) temperature. 
In Tab. IHIl we report the corresponding e- and p- values of the two structures identified by the 
evolutionary algorithm. None of the additional crystals obtained from the MC simulations, 
neither the FCC nor the Ts lattices, are identified as global minima of the Gibbs free energy 
at zero temperature. Nonetheless, both structures, at high and low pressure respectively, are 
among the best configurations found by the optimization technique: they are both identified 
as low- lying local minimum of the Gibbs free energy, in their corresponding range of pressure. 

To conclude, we also consider the regular tetrahedral case corresponding to g = 109.47°. 
The phase diagram of such a system is reported in Ref.— ; it shows the presence of three 
different FCC-like phases: a low temperature face-centered non-cubic structure, an FCC 
phase, and a plastic FCC crystal (FCC P ). The EA approach only identifies the first one as 
the global enthalpy minimum at zero temperature; however, the FCC lattice is identified 
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as a local minimum of the enthalpy landscape^. In contrast, the MC approach is able to 
identify the latter two FCC-like phases. In Fig. HI we compare the frequencies of appearance 
of the packed structures found via the variable box shape simulation technique over a large 
temperature range. Among the packed structures, we also identify a BCC like phase, which 
has a wide region of stability in the phase diagram^. Fig. H] shows that the variable box 
shape simulation technique is able to properly take into account the effect of temperature; 
indeed at the highest temperatures investigated, FCC P is observed with increasing frequency. 
On the other hand, with the MC approach, we do not observe the non-cubic structure, whose 
region of stability is confined to very low temperatures^. 

V. CONCLUSIONS 

In this paper we have employed the MC NPT variable box shape simulations to predict 
candidate structures for several patchy particle models proposed in the literature. We have 
determined the structures for two patchy hard-sphere models, i.e. the Kern-Frenkel model 
with four patches and the "sticky spots" model with six patches. For the Kern-Frenkel model, 
we find all the stable phases as previously predicted in Ref.— thereby giving confidence 
in the MC method. Moreover, we find a BCT and two HCP phases which stability should 
be determined by free energy calculations. For the sticky spots model, we have successfully 
predicted several candidate structures. To the best of our knowledge, the crystal phases of 
this system have not been studied before. 

In addition, we have compared crystal structures predicted by MC NPT variable box 
shape simulations with results obtained via an evolutionary algorithm approach for various 
patchy particle systems interacting via continuous pair potentials. From our findings, it 
appears that neither method is significantly better than the other, and that the most ap- 
propriate method for a given system depends on the characteristic features of the problem. 
Concerning the relative efficiency in finding the solid structures, approximately 90% of the 
EA runs converge to the same minimum of the enthalpy landscape, irrespective whether it 
corresponds to an open or a packed configuration. On the other hand, MC runs show two 
different percentage values in the two cases: almost 100% in the high/intermediate pressure 
regime and around 75% in the low pressure regime. Moreover, the computational costs of 
the two approaches differ by less than one order of magnitude on a real time scale. In the 
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following, we briefly highlight a few of our considerations on the comparison between the 
two numerical approaches: 

(i) For low temperatures (T ~ 0), the two techniques produce virtually equivalent crystal 
structures, when a selection of the lowest-lying local minima on the enthalpy land- 
scape identified by the evolutionary algorithm is taken into account: structures listed 
among the energetically most favorable ones according to the evolutionary algorithm 
technique, despite not being a global minimum at zero temperature, can be possibly 
thermodynamically stable at finite temperatures^. This encouraging fact is a hint on 
the reliability of both methods. Structures that are identified as T = stable phases 
can be then considered as good candidates at finite (even though low) temperature. 

(ii) Further, the MC method can predict candidate crystal phases at finite temperature, 
unlike the evolutionary algorithm, which is bound to T = imposed by conceptual 
and computational limitations. Indeed, on increasing the temperature beyond the 
chosen range (T > 0.2), we identify FCC plastic crystal phases for all the discussed 
models. Moreover for the regular tetrahedral Lennard- Jones patchy system we are 
able to identify a FCC crystal structure which is stable only at finite temperature 
and cannot be found via the EA approach. On the other hand, at the temperatures 
relevant for patchy systems the simulations are more likely to get kinetically trapped 
in gel-like states or non-competitive local minimum configurations, which have to be 
discarded. This problem is much easier to handle within the evolutionary algorithm 
by using suitably designed "population control" operations". 

(iii) For systems with discontinuous interaction potentials, the MC method has the advan- 
tage that it can be applied directly, while for use with an evolutionary algorithm either 
suitable approximations (by smoothening the potential) or cumbersome methodologi- 
cal implementations^ are needed. 

We stress again, that the thermodynamic stability of the crystal structures predicted 
by both methods is not guaranteed and has to be checked by full free-energy calcula- 

As the two methods covered here have both advantages as well as shortcomings, we list 
potential improvements in the following. A possibility to overcome kinetic trapping in gel- 
like states for variable box-shape simulations is to combine this method with moves that 
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correspond to larger leaps in configuration space, comparable to mating and mutation op- 
erations of an evolutionary algorithm; this would move the method further away from being 
a thermodynamic approach into the direction of an optimization technique. The evolution- 
ary algorithm, on the other hand, could be augmented with especially effective mutation 
steps, that are based on short MC runs, as suggested in Ref.— . Ultimately, even a hybrid 
approach, incorporating the advantages of both methods, is conceivable. Another desirable 
amendment to EAs would be using free energy calculations based on lattice dynamics^ 3 - in 
order to estimate the competitiveness of candidate structures at finite temperature already 
during the run of the algorithm; it has to be noted though, that such an approach demands 
the interaction potential to meet even stronger criteria (continuous second derivatives) and 
is conceptually rather involved and computationally expensive. 
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structure 


e 


P 


/ 


HCP (a) 


-4/3 


1.37 


8% 


HCP (b) 


-4/3 


1.37 


4% 


FCC (c) 


-2 


1.18 


53% 


BCT (d) 


-2 


1.25 


26% 


DC (e) 


-2 


0.6 


30% 


DH (f) 


-2 


0.6 


23% 



TABLE I: Overview of the crystal structures shown in Fig. [T]for the four patch Kern-Frenkel model. 
For each structure, the corresponding values of the average energy per particle e (in units of e), 
the average density p (in units of cr~ 3 ), and the frequency of appearance / (%) are reported. The 
upper part of the table contains the four structures found for high/intermediate pressure values 
(over a total of 160 different simulations started in parallel with different initial conditions); the 
part of the table below the horizontal line contains the two open structures found at low pressure 
values (over a total of 70 different simulations). The densities of the DC and the DH crystals are 
reported with a one-digit precision (vs a two-digits precision in the case of close-packed crystals) 
since, for a square well attraction of range 5 = 0.24cr, the density of an open bonded structure can 
vary up to 10% without causing an additional cost in energy. 
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structure 


e 


P 


/ 


HCP (a) 


-2.25 


1.37 


1% 


FCC (b) 


-3 


1.37 


3% 


FCC (c) 


-2.5 


1.37 


7% 


FCC (d) 


-2 


1.18 


38% 


BCT (e) 


-3 


1.33 


8% 


SC (f) 


-3 


0.87 


44% 


Hcl (g) 


-3 


0.80 


7% 



TABLE II: Overview of the crystal structures shown in Fig. [2] for the sticky spots model with six 
patches per particle. For each structure, the corresponding values of the average energy per particle 
e (in units of e), the average density p (in units of c -3 ), and the frequency of appearance / (%) 
are reported. The upper part of the table contains the five structures found for high/intermediate 
pressure values (over a total of 315 simulations); the lower part of the table below the horizontal 
line contains the two open structures found at low pressure values (over a total of 225 simulations). 
The densities of the open structures are reported with a two-digits precision (as in the case of 
close-packed crystals), since the change in density - without additional energy costs - of an open 
bonded structure can amount to 6% at maximum in the case of a square well attraction of range 
5 = 0.119a. 
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Geometry 


structure 


MC 


EA 
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e 
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P 


minimum 


I 


HCP (a) 


-0.83 


1.23 


17% 


-1.19 


1.34 


9.1 


LMi 


I 


FCC (b) 


-1.03 


1.18 


80% 


-1.23 


1.32 


9.1 


GM 


I 


HxDl-II 








-1.18 


1.35 


6.0 


GM 


I 


HxDl (c) 


-1.17 


0.91 


24% 


-1.68 


0.98 


1.1 


LM 5 


I 


HxDl-I 








-1.84 


1.04 


1.1 


GM 


I 


HcDl (d) 


-1.50 


0.57 


50% 


-2.00 


0.65 


0.1 


GM 


II 


HCP (e) 


-1.13 


1.30 


33% 


-1.69 


1.33 


6.0 


GM 


II 


FCC (f) 


-0.88 


1.18 


66% 


-1.66 


1.32 


6.0 


LMi 


II 


HxSDl (g) 


-1.63 


1.02 


45% 


-2.00 


1.11 


0.1 


GM 


II 


Ts (h) 


-1.54 


0.88 


30% 


-2.00 


1.03 


0.1 


LM 2 



TABLE III: Overview of the crystal structures for the orientational Lennard-Jones model with 
four patches per particle. Structures identified by both the Monte Carlo (MC) and the evolutionary 
algorithm (EA) approach are listed according to the labels of Fig. [2 for the structures that are 
obtained only from the EA method we refer to Ref.— . The first column of the table specifies 
the patch geometry: case I corresponds to g = 90°, while case II corresponds to g = 150°. The 
values of e (in units of e), p (in units of cr~ 3 ) and / for the candidate structures obtained via the 
MC simulations are shown in the third, fourth and fifth columns, respectively. The sixth and the 
seventh columns report the e- and p- values of the structures as identified by the EA approach; these 
values are obtained at the pressure reported in the 8th column. When a given structure is (locally) 
optimized for a different pressure value, e and p can change by about 5%. The ninth column gives 
information about the location of the corresponding minimum on the Gibbs free energy landscape 
at T = 0: GM corresponds to the global minimum and LM n corresponds to the nth local minimum 
after the global one. For both patch geometries, the upper part of the table includes structures at 
high/intermediate pressure values (41 MC simulations for case I, 15 MC simulations for case II) 
and the lower part contains open structures at low pressures (50 MC simulations in case I, 20 MC 
simulations in case II). 
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FIG. 1: (Colors online) Representative parts of the crystal structures identified for a particular 
realization of the Kern-Frenkel model with four patches per particle^. The colors of the particles 
represent the number of bonds per particle: light red = four bonds (i.e. fully bonded), magenta = 
three bonds, dark red = two bonds. Sticky sites are colored in green. Packed structures: (a) and 
(b) HCP crystals with same average energy per particle and same average number density, (c) fully 
bonded FCC crystal, (d) fully bonded BCT crystal. Open structures: (e) part of a DC crystal and 
(f) part of a DH crystal. The yellow circles highlight the different intra-layer bonding in the DC 
and the DH structures (in order to highlight the patch positions, a different particle/patch size 
ratio is chosen; in each circle two particles, one on the top of the other, are reproduced): pairs of 
particles forming bonds between six-fold rings have a staggered conformation in the DC crystal, 
and an eclipsed conformation in the DH crystal. Tab. [J contains the e- and p- values, as well as the 
/-values, of the respective candidate structures. 
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(e) CO (g) 



FIG. 2: (Colors online) Representative parts of crystal structures identified for the sticky spots 
model with six patches per particle^. The colors of the particles represent the number of bonds 
per particle: light red = six bonds (i.e. fully bonded), blue = five bonds, turquoise = four bonds. 
Sticky sites are colored in green. Packed structures: (a) HCP structure, (b) fully bonded FCC 
crystal, (c) and (d) partially bonded FCC crystal distinguished by different e- and p-values, and 
(e) fully bonded BCT. Fully bonded, open structures: (f) SC crystal and (g) honeycomb layers 
(Hcl) (planes of six- fold planar rings). Tab. HT1 contains the e- and p- values, as well as the /-values, 
of the respective candidate structures. 
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FIG. 3: (Colors online) Representative parts of crystal structures identified for the orientational 
Lennard-Jones model with four patches per particle introduced in Ref.— . The selected patch ge- 
ometries are displayed in the rectangular box on the left hand side: case I corresponds to the 
compressed tetrahedron with g = 90°, and case II corresponds to the elongated tetrahedron with 
g = 150°. Representative parts of the observed structures for case I : (a) HCP crystal, (b) FCC 
crystal, (c) hexagonal double layers (HxDl), and (d) honeycomb double layers (HcDl). Represen- 
tative parts of the observed structures for case II: (e) HCP crystal, (f) FCC crystal, (g) hexagonal 
staggered double layers (HxSDl), and (h) tetrahedrally arranged lattice (Ts). Tab. Ill II contains the 
e- and p-values, as well as the /-values, of the respective candidate structures. 
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FIG. 4: Overview of the packed crystal structures for the orientational Lennard-Jones model with 
four patches per particle arranged in a regular tetrahedral geometry {g = 109.47°). The histogram 
reports the /-values (normalized to unity) of each identified structure at different temperatures (in 
units of e/ks) indicated at the bottom of each bar. 
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